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Abstract 

We derive an analytic expression for the scalar one-loop pentagon and hexagon functions 
which is convenient for subsequent numerical integration. These functions are of relevance in 
the computation of next-to-leading order radiative corrections to multi-particle cross sections. 
The hexagon integral is represented in terms of n-dimensional triangle functions and (n-|-2)- 
dimensional box functions. If infrared poles are present this representation naturally splits into 
a finite and a pole part. For a fast numerical integration of the finite part we propose simple 
one- and two-dimensional integral representations. We set up an iterative numerical integration 
method to calculate these integrals directly in an efficient way. The method is illustrated by 
explicit results for pentagon and hexagon functions with some generic physical kinematics. 



1 Introduction 



The growing potential of current and future high energy cohiders permits an increasingly detailed 
and precise study of particle physics and will push the energy frontier to the TeV scale before the 
end of this decade. Experiments at the Fermilab Tevatron, the LHC at CERN and the planned 
e+e" linear collider will produce vast amounts of data. The prospect of unprecedented experimental 
statistics calls for a corresponding level of precision on the theoretical side which generally requires 
higher-order predictions for the underlying multi-particle interactions. To achieve this objective is 
particularly important for the discovery of physics beyond the Standard Model, where a precise 
description of Standard Model backgrounds is crucial and necessitates next-to-leading order (NLO) 
calculations for a variety of multi-particle search channels. The resulting reduction of unphysical 
scale uncertainties allows for precise quantitative predictions, thus significantly improving the rather 
qualitative leading order descriptions. For background predictions the importance of higher-order 
effects is further enhanced by severe experimental cuts, because the latter typically select peripheral 
phase space regions, which are particularly sensitive to NLO effects. 

While for partonic and electroweak 2^2 processes advanced theoretical knowledge about NLO 
calculations exists, our understanding of processes with more final-state particles has not yet reached 
maturity. In recent years, many relevant 2—^3 amplitudes and processes involving a small number 
of scales (as typical in QCD) have been calculated at NLO. However, results for multi-scale 2 — 3 
processes are rare, and the step to 2 ^ 4 amplitudes again entails a significant increase in complexity. 
Currently, 2 4 results are only known for certain supersymmetric gauge theory amplitudes and 
the Yukawa model Q . 

NLO calculations for a 2 ^ process are typically organized by evaluating the real emission 
2 ^ + 1 part and the virtual corrections to the 2 N process separately]^. In the latter, 
divergences are regulated by working m n — A — 2e dimensions. The analytic evaluation of the 
virtual part requires the reduction of tensor integrals to scalar integrals, which is well understood 
for arbitrary 2 — > TV processes [Q, ||, ^ . Reduction formulas for scalar integrals which relate general 
finite iV-point functions to box integrals have been formulated in for the infrared finite case. For 
massless integrals with N < 7, similar reduction formulas have been derived in n dimensions in 
The generalization to arbitrary N has been considered in |9| for massive integrals and in |^ for 
the massless case. 

As will be pointed out below, the reduction approach leads to a natural separation between IR 
divergent and finite contributions. The IR divergences can be collected in triangle functions with 
massless propagators which have simple analytic representations. The finite remainder of the scalar 
integrals can be expressed fully analytically in terms of dilogarithms (and simpler functions) |l0|| . 
However, the large number of kinematic invariants leads to a huge number of dilogarithms, many 
thousands in the case of the massive hexagon function. A numerical evaluation at this level typically 
leads to large cancellations in certain kinematic regions and thus to numerical instabilities Jill . Hence 
a numerical evaluation at an earlier stage may be equally good for practical purposes, if not better, 
and this is what we suggest in this paper. 

An alternative numerical approach to the evaluation of Feynman diagrams which deals with 
the computation of multi-leg one- loop diagrams has been presented in It is based on 

the Bernstein-Tkachov theorem |Q. The basic idea is to rise the power of negative exponents 
of kinematic functions by using the Bernstein-Tkachov relation. This leads in principle to better 
behaved integrands. However, an explicit result for the hexagon function has not been given yet. 

In this paper we derive a representation of the hexagon function in terms of 20 n-dimensional 
triangle functions and 15 (n-|-2)-dimensional box functions. For completeness we also provide the 
pentagon function in terms of 10 triangle and 6 box functions. By using sector decomposition 
and one explicit integration we derive a one-dimensional integral representation for the triangle 

^An exception is the purely numerical approach of Soper which combines real and virtual corrections on the 
integrand level. 
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and a two-dimensional parameter integral for the box. We study the analytic structure of these 
representations in some detail. The virtue of these parameter integral representations lies in the fact 
that the singularity structure is quite transparent. Hence, as a by-product, we derive the well-known 
fact that for physical kinematics only integrable singularities of logarithmic and square-root type 
are present at one loop. 

Further, we describe a numerical integrator that facilitates a fast and accurate evaluation of 
these "atoms" of our representation. Finally we give numerical examples for our approach. We only 
consider the case where all propagators are massive here, but also comment on how to generalize 
our approach in the presence of IR divergences. 

The paper is organized as follows. In Section 2 we derive representations for the n-dimensional 
box, pentagon and hexagon functions in terms of rt-dimensional triangle and (n-|-2)-dimensional box 
functions. In Section 3, one- and two-dimensional integral representations are given for the triangles 
and boxes, together with a detailed discussion of the singularity structure of the integrands. In 
Section 4 we outline the numerical evaluation of these integral representations. Examples for our 
procedure are provided in Section 5. The article closes with a summary and an outlook. 



m 



2 Reduction of pentagon and hexagon integrals 

In this section we will provide explicit expressions for the n-dimensional hexagon, pentagon and 
box functions in terms of n-dimensional triangle and (n-t-2)-dimensional box functions. If infrared 
divergences are present, which manifest themselves in terms of poles in 1/e = 2/(4 — n), this choice 
of building blocks naturally splits these functions into infrared divergent and finite pieces. After 
separating the IR divergent triangles from the representation, only integrals which can be evaluated 
in n = 4 dimensions remain. To fix our conventions we define the n-dimensional A^-point function 
as 

/ \ f d k 

I]y{pi, . . . ,PN,mi, . . . ,mNj = 

The kinematic information is contained in the matrix S which is related to the Gram matrix G by 

Ski = -{n-rkY +mf +ml^{Gki-vi-Vk) (2) 
Gki = 2rk-ri, Vk = Gkk/2 - ml , k,l^l,...,N 

The vectors Vj = pi + . . . + pj are sums of external momenta. In physical applications the external 
momenta span the 4-dimensional Minkowski space. We consider here only the case where any four 
of the A'' external vectors are linearly independent. Further we assume that the kinematics is such 
that the anomalous threshold, i.e. the leading singularity of the A^-point function, is not probed. 
We recall that the leading singularity of a Feynman integral in parameter space corresponds to a 
vanishing denominator of the parameter integral while all values of the Feynman parameters are 
nonzero p^ . This means that the matrix S has to have a zero eigenvalue or simply that det(S') = 0. 
The generalization of our approach to this exceptional case is briefly discussed below, but let us 
first focus on non-exceptional kinematic configurations. In this case every A^-point function with 
Af > 6 can be expressed in terms of pentagon integrals Q . The pentagon functions, in turn, can be 
expressed in terms of box functions up to a term which vanishes in the limit n — > 4. For the purpose 
of this paper, we only give the reduction formula for the case A^ < 6. It is derived in the same way 
as in |5| , where a simple derivation for the massless case and general A^ is explicitly given. 

In = flBkr^,_,,k + {N ^n-l)^^I-+^ , det(^) ^ 0, (3) 
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N 



(4) 



1=1 



Hence the A^-point function decays into a sum of (A^-l)-point functions, which are obtained by 
pinching propagators 



rn 



,;_n/2 



(5) 



and a remainder term which is proportional to the Gram determinant, det(G), times the (n+2)- 
dimensional A'^-point function. The latter turns out to be infrared finite even in the massless case, 
as can be seen by power counting. The reduction coefficients Bk are defined through the inverse of 
the kinematic matrix S. Note that the formula is valid for general dimension, as long as the external 
momenta are non-exceptional. The rank of the Gram matrix is min(4, A'^ — 1), whereas the rank of 
the matrix S is min(6, iV), such that the inverse of S does not exist for > 6. However, we stress 
that reduction formulas do not rely on the regularity of the matrix S, since in the case of a singular 
S one can follow the lines of |^ , where the concept of a pseudo- inverse to a matrix was used to deal 
with the case A^ > 6. The case A^ < 6 with det(5) — can be treated analogously. 

In the generic case of non-exceptional momenta the hexagon function decays into six pentagons 
without rest 



Iq{si2,S23, S34, S45, S56, Sgl, S123, S234, S345J Sli ■ ■ ■ iSq, "li, . . . , mg) — 

Bi /^(si23, S34, S45, S56, S345, S12, S3, S4, S5, se,m2, m3,mi, TO5, uiq) 

+B2 /?(S234, S45, S56, S61, S123, S23, S4, S5, Sg, Si, ^3, 7714, TO5, UIq, TOl) 
+^3 -^5'(S345, S56, Sgi, S12, S234, S34, S5, Sg, Si, S2, m^, ms, TOg, mi, m2) 

+B4 /^(si23, S61, S12, S23, S345, S45, S6, Si, S2, S3, ms, me, wi, m2, m3) 

+B5 /"(S234, S12, S23, S34, S123, S56, Si, S2, S3, S4, TOg, mi, TO2, m3, TO4) 
+i?6 -^^(S345, S23, S34, S45, S234, Sgi, S2, S3, S4, S5, TOi, m2, TO3, m4, TO5) 



(6) 



Bi 



(6)n 
kl ) 



1=1 



The kinematic matrix is 



Sif+ml + mf 



0(6) 



5(6) 



The kinematic invariants are defined as s,- 



(fc,/ = l,...,6) 
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S234 
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S3 





S4 


S45 


S123 
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S45 


S5 





S6 
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S123 
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Pj, Sij, 



iPi + Pj 



Similarly, the general pentagon integral can be written as 

-fr(si2, S23, S34, S45, S51, si, . . . , S5, mi, . . . , TO5) = 

Bi (S45, S34, S12, S3, S4, S5, TO2, m3, TO4, ms) 

+B2 /4'(s5i, S45, S23, S4, S5, Si, TO3, m4, m5, mi) 
+B3 /4 (si2, S51, S34, S5, Si, S2, TO4, ms, mi, m2) 
+B4 /4 (s23, S12, S45, Si, S2, S3, ms, mi, 1712, ms) 
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+^5 -fr(s34, S23, S51, S2, S3, S4, mi, 1712, m^, mi) + 0{e) 



(7) 



here Bk 



E 

1=1 



(5) 



kl 



S 



(5) 
kl 



(fc,; = i,...,5) 



^(5) = - 



/ S2 

S2 

S23 S3 

S51 S34 S4 

V Sl S12 S45 S5 



S23 S51 Sl \ 

S3 S34 S12 

S4 S45 

S5 

y 



In principle, if no infrared divergences are present, there is no need to reduce any further, as analytic 
formulas for the finite 4-dimcnsional box integral exist. The boxes can be expressed in terms of a 
large number of dilogarithms ||l6|. However, as is well known these representations are not 
unproblematic from a numerical point of view, since large cancellations between the dilogarithms 
occur in certain kinematic regimes. For this reason, and also in view of the general case which 
includes infrared divergences, it will turn out to be useful to reduce the box integrals further. This 
allows for a natural separation of the IR-singular and finite terms. 
The reduction formula for boxes reads 



^r(si2, S23, Sl,..., S4, mi, . . . , 7714) = 

Bi ^3 (si2, S3, S4, m2,m3, nii) 

+ B2 13(323, 34, Sl, 7713, 1714, mi) 

+B3 13(312,31, 32, m^, mi, ma) 
+i34 /J (S23 , S2 , S3 , mi , m2 , m3 ) 

+ (n - 3)(Bi +B2 + B3 + Bi) /r+^(si2, S23, Sl, . . . , S4, mi, ... , m4) (8) 



Bk 



5(4) 



1=1 

/ S2 

S2 

S23 S3 

\ Sl Sl2 



S. 



(4) 
kl 



c.(4) 
'^kl 



+ ml 



(fc,/ = l,...,4) 



S23 
S3 



S4 



Sl \ 

Sl2 

S4 

/ 



Note that det(G('*)) = -{Bi + B2 + B3 + B4) det(5'('^)). 

By applying the reduction formulas iteratively one finds a representation of the pentagon integral 
in terms of 10 triangle and 5 (n+2)-dimensional box integrals. It can be written in the following 
compact, cyclically symmetric form 



4" 



Bi5 



Tn+2 

4,1 



B25)in' 



(B1B12 + B2B2l)l3 12 + (Bli3l3 + -63-631)13 J3 + Bi(Bi2 + -B13 + Bii 
-{B2B23 + B3B32)l3 23 + (B2B24 + -B4-B42)-?3,24 + -62(-S21 + -B23 + -B24 
-{B3B34 + 54543)73 34 + (B3B35 + 55553)73 35 + 53(531 + 532 + 534 + 535)74 
-(B4B45 + 55554)73 45 + (B4B41 + BiBi4)I^i^ + 54(541 + 542 + 543 + 545)74 
(55551 + 51515)73 ^5 + (B5B52 + 52525)73 25 + 55(551 + 552 + 553 + 



n+2 
4 



r> \ Tn+2 
554j74 5 



(9) 



Here the 5^ are the reduction cocfhcients of the pentagon integral (Q) and 5jj is the jih reduction 
coefficient of that box integral which stems from the ith pinch of the pentagon integral. Note that 
Bij 7^ Bji. On the other hand, the triangles, which result from double pinches of the pentagons, are 
symmetric: 73 ,y — I3 ji. Analogously, one finds for the hexagon integral a representation in terms 
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of 20 triangle and 15 (ri+2)-dinicnsional box integrals: 
^6 — 

|[Bl(i?12i3l23 + -813^132) + B2{B2lBi23 + B23B231) + B^{B^iBij,2 + -B32-S231)] ^3*423 + ^ C.p.j 

+ \[Bi{Bi2Bi2i + B14S142) + -B2(-B2ii?2i4 + S24-B241) + Bi{BiiBii2 + -642^421)] ^3j24 + 5 c.p. I 

+ 1 [-^1 (^13-^134 + i?i4Si43) + i33(i?3ii?3i4 + ^34-6341) + 54(541-6413 + -643^431)] ^3^134 + 5 C.p.j 
+ (^13^135 + ^15^153) + i?3(i?3ii?3i5 + B^<^B^^i) + B<^{B<^iB<^i^ + -653^531)] -^3^135 + 1 C.p.j 

+ [{BiBi2 + B2B2l){Bi23 + Bi24 + Bi25 + Bi2fi)l2X2 + 5 C. p. j 
+ {(BlBi3 + B3S3l)(Bi32 + B134 + B135 + Si36)/!l3 + 5 C.p.j 
+ |(i?l_Bi4 + i?4i?4i)(_Bi42 + i?143 + ^145 + Bi4q)I2j^I + 2 C.p.j 

(10) 

Here the Bi are the reduction coefficients of the hexagon integral (^. The Bijk = Bjik are a 
shorthand for the reduction coefficient of the fcth pinch of the box integral I^^j and c.p. means 
cyclic permutation of the indices of the _B's and the indices which define the pinches. The prob- 
lem of calculating the pentagon and hexagon integrals is now reduced to the calculation of lower 
point integrals and reduction coefficients. As pointed out above, a complete analytic expression is 
possible but not necessarily of advantage. Such an analytic expression contains a huge number of 
dilogarithms, and nontrivial numerical cancellations occur during evaluation. Since one has to rely 
on numerical integration at some stage of the calculation of most physical cross sections anyhow, 
a direct numerical evaluation of the scalar integrals in parameter form is more than adequate for 
practical applications. In our approach, all one has to do is to provide stable and sufficiently fast 
numerical integrators for the finite 4-dimensional triangle and 6-dimensional box integrals. 

Numerical instabilities typically arise from terms with opposite signs and denominators that 
approach zero. The denominators that occur in our reduction are the determinants of the kinematic 
matrices S from the different reduced hexagon, pentagon and box integrals. Thus, the critical 
points are the normal and anomalous thresholds fl^ of the corresponding scalar graph. Near these 
thresholds one finds that the reduction coefficients fulfill (approximately) additional constraints. 
This can be exploited to achieve stable groupings of terms in the respective limits. 

3 Integral representations of triangle and box functions 

In this section we first derive integral representations of the triangle and box functions which are 
appropriate for direct numerical integration. Then we thoroughly analyse the singularity structure 
of the integrands. 

Using standard Feynman parametrisation, the parameter representations of the n-dimcnsional 
3- and 4-point functions are 

jn/ 2 2 2\ _ [ "^"^ 1 

i3isi,S2,S3,mi,m2,m3j - j ^^n/2 ^^^k - r,f - ml][{k - r2Y - mlW - 

00 

= -r(3-n/2) / d^x^(l^xi23) /''^'^Xj" (11) 
Ttti = (-si)xia;3 + (-S2)a;ia;2 + (-S3)a;2a;3 + a;i23(a;imi + a;2m2 + a;3m3) - i(5 

2;i23 = Xi + X2 + 
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and 



I2{S12, S23, Si, S2, S3, S4, m^, mj, TO3, 7714) 



z7r"/2 [(k — ri)2 — TO^][(fc — — 7712] [(/c — ra)^ — TO3][fc2 — TO4] 

00 

r(4-n/2) / d^xJ(l-xi234) /j^^^^r_7/2 (12) 

i (^Box) 



Tbox = (-Si2)x2a;4 + (-323)0:1X3 + [-si)xiXi + [-S2)xiX2 + (-S3)x2a;3 + (-S4)x3a;4 

+Xi234:{ximl + X2ml + X3ml + X4^ml) (13) 

Note that we will need the case n — 6 — 2e for the box in the following. 

In special cases, e.g. if the integral has a high symmetry due to equal scales or vanishing 
kinematic invariants, a non-symmetric transformation of Feynman parameters typically leads to 
simplifications. However, here we want to deal with general situations, so no Feynman parameters 
should be preferred and a symmetric treatment of the problem seems adequate. This is achieved 
by splitting the A^-dimensional {N — 3,4) parameter integrals into N sectors by the following 
decomposition 

1 = 0{xi > X2, . . .^Xn) + 0{X2 > Xi,X3, . . .,Xn) + . . . + 9{XN > Xi,.. .,XN-i) (14) 

This approach will turn out to be useful for numerical integration later. In addition, it is very 
convenient in the presence of IR divergences |]l^ . The box integral then decays into a sum of four 3- 
dimensional parameter integrals, while the triangle integral decays into a sum of three 2-dimcnsional 
parameter integrals. Focusing on the last sector as an explicit example, one finds after the variable 
transformation tj — Xj/x^ {j = I, . . . , N — 1) 

1 

c" / 2 2 2\ f j+ (1 + ^1+^2) ~" 

Jt«(S1, S2, S3, TOi, 7712, 7713) = / dtidt2—— .^/2 ^^^1 

( ^Tri ] 

TTr^ = (-Sl)tl + (-S2)tl^2 + (-33)^2 + (1 + tl + t2)(tl™? + ^2^2 + ml) (16) 

1 

'S'box (■512, S23, Si, S2, S3, 34,777,1,7772,7773,777.4) = dtidt2dt3 l3-n% (1'^) 



Fbox = (-312)^2 + (-S23)il^3 + (-Sl)<l + (-S2)tlt2 + (-33)^2^3 + (-34)^3 

+ (1 + tl + t2 + H){t\m\ + t2m\ + ^3771^ + mX) (18) 

Since in the Euclidean region both functions Tbox and Ttti are strictly positive, bounded functions, 
it is clear that the respective integrals can easily be computed numerically. We will see shortly 
that the above representations allow for a transparent discussion of the singularity structure of the 
integrals for physical kinematics. The box and triangle functions are given in terms of the sector 
functions ( IT^ and ( |l7|) as 

^?(si,S2,S3,777?,77i^,m2) = -r(3 - 77/2) X 

on I 2 2 2\ , on r 2 2 2\ , on l 2 2 2\ 

'3t«(S2, S3, Si, 7772, 7773, "^l) + ^Tri{S3,Si, S2,m^, ^i, 77I2) + Srp„[Si, S2, S3, 77li, 7772, "^3) 

(19) 
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Tn+2 



(si2, S23, Si, S2, S3, S4, mf, ml, ml, ml) = r(3 - n/2) x 



01+2/ 2 2 2 2\ I cn+2/ 2 2 2 2\ 

^Box (■523, S12, S2, S3, Si, Si, TO2, TTig, TO4, mj) + S (si2, S23, S3, S4, Si, S2, m^,m^, m^,m2) 

I c'ra+2/' 2 2 2 2\ , c'n+2,- 2 2 2 2n 

+ !^Box i*23, Sl2, S4, Si, S2, S3, 7714, TOi, mj, mg) + (Sl2 , S23 , Si , S2 , S3, S4, ^i, mj, Wg, 7714) 

(20) 

If no infrared divergences are present one can set 77 = 4 in both formulas. Explicit formulas for IR- 
finite integrals in 77 = 4 dimensions are known for triangles and boxes and can be found in pq, 
For the case of vanishing internal masses a list of box and triangle integrals may be found in 
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PI- 
I- 

To the best of our knowledge no complete list of all mixed cases is given in the literature. We note in 
this respect that our representation of 4-dimensional box integrals in terms of triangle functions and 
6-dimensional boxes is a good starting point, as the infrared divergences are exclusively contained 
in triangle functions in this case. 

In both integrals, ( |l5| ) and ([l7|), the exponent of the kinematic functions is —1. Note that 
high negative powers of kinematic functions are typically the reason for the failure of a stable 
numerical evaluation of multi-leg or -loop Feynman diagrams in the physical region. In this respect 
our representation is well suited for a numerical approach. As the kinematic functions are quadratic 
in each parameter, one integration can be done analytically. For the triangle one gets 



1 

'S'th'*(s1, S2, S3, 777^,7772, 7773) = / dtidt2 , ^ , , 

J U + ii 



t2) Atl + Bt2 + C -i5 



dU 



2A 



with 



"log(2A + B-VR)- log(5 -VR)- log(2A + B + T) + logjB + T) 

T+y/R 

logi2A + B + VR)- log(B + Vi!) - log(2A + B + T)+ log(g + T) 

t-Vr 



A = ml 

B = {ml + ml — S2)ti + ml + ml — S3 

C ~ mftf + {mf + ml — Si)ti + ml 

R = B^ - AAC + i5 

T = 2A{l + ti)-B 



(21) 



(22) 



For the 6-dimensional box integral one finds similarly 

1 



^B:.^(s12, S23, Si, S2, S3, S4, 777^ 7^^ 7^^ ml) ^ j dtMH ^ ^ ^ ^ g _ 



/■ , , 4^2 f AA^fR 



1 



1 



(T - Vi?) 



\\{2A + B + T){B + T){T^ - R) 
log(2A + B- VR) ~ log{B - VR) - log(2A + B + T) + log(B + T) 

log(2A + B + Ve) - log(B + VR)- \og{2A + B + T)+ log(B + T) 



(23) 



7 



witli0 



B = {ml + ml — 523)^1 + + "^3 ^ ■53)^2 + ml + ml ~ 34 

C = (-512)^2 + (-si)ti + (-52)^1^2 + (1 + ii +<2)(m| + m^ii +TO^i2) 

R = - 4:AC + iS 

T = 2A{1 + h + t2) - B 



(24) 



The critical points for the numerical evaluation of both integrals are vanishing denominators and 
logarithms with arguments tending to zero. Before analysing the singularity structure of the inte- 
grands we explicitly separate imaginary and real part. It is useful to distinguish the cases R > 
and i? < in this respect. If i? < 0, then C > and A + B + C > 0, no imaginary part is present. 
We find 



Sxri (si,S2,S3,mi,m2,m3) 



dti 



4A 



T^-R 



{ [log(2A + B + T)- log(B + T)] + e{R < 0) 



log(C) - \og{A + B + C) 



T 



-e{R > 0) 



: arctan 



-R 



B 



arctan 



2A + B 



n e{B <{)<2A + B) 



T-^/R 
2^/R 



log {\2A + B- v^l) - log (\B - Vr\) - inOiB <Vr<2A + B) 



2y/R 



log (\2A + B + V^l) - log (|B + V^l) + iTr0{B < -V^ <2A + B)) } 



(25) 



and 



^Boxisi2, S23, Si, S2, S3, Si, ml, ml, ml, ml) 



dtidt2- 



16^2 



A{T^ - R) 



{2A + B + T){B + T) 

-.(i^<o)r^^°s(^)-^°^^^ + ^ + ^) 

T^ + R 



(r2 - i?)2 

u 

+ T [\og{2A + B + T)- \og{B + T)] 



-R V 



^e{R > 0) 



( arctan 



^-R 



B 



-R 



{T - VRf 
4%/i? 



arctan | , '^ j + tt e{B < Q < 2A + B) 
{log {\2A + B- \/^|) - log (\B - Vr\^ -iTT 0{B <Vr<2A + B) 



(T+VR)^ 
4 



(log (\2A + B + VR\^ - log (\B + + in e{B < -Vr < 2A + B))] } 



(26) 



We define the step function to be 1 if its argument is true, and else. Let us now investigate 
the singularity structure of the integrands. At first sight, dangerous, possibly singular denominators 

^We do not introduce new symbols A, B, C, R, T for the box. Which definition applies is always clear from the 
context. 
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are present. First note that the hmit T — > with simuhancously i? — > is unproblematic since 
B > 2A > and C > A > in this case. The integrands then are bounded and positive definite, 
as can be seen by looking at the starting expressions in (|2^) and (p3|). Further, a short calculation 
shows that the limits T ±\/i? with R > are finite, so the integrand is also non-singular in this 
limit. The last denominator which can vanish is \/±i?. In the limit ±R the integrands in ( p5| ) 
and ( p6| ) behave as 

^n6{B<Q<2A + B) l^AE^S- + ^ + finite (27) 

L \ —R \ R J 



As i? is a quadratic form in the integration variables, we have an integrable singularity of square- 
root type. The integrand also exhibits logarithmic singular behaviour whenever an argument of 
a logarithm goes to zero, i.e. at the boundaries of the regions where the integrand develops an 
imaginary part. Hence, it is necessary to have i? > in order to produce a logarithmic singularity. 
Three regions which lead to an imaginary part can be distinguished: 

Region I: A + B + C > 0, -2A <S<0,C>0^(B< ±\/R <2A + B). 

Region II: A + B + C > 0, C < ^ (B < \/^ < 2yl + B) and not (B < -^/R <2A + B). 

Region III: A + B + C<0,C>0^{B< -y/R <2A + B) and not {B <y/R<2A + B). 

Region I is an overlap region where the imaginary part has two contributions. In regions II and III 
only one of the ^-functions in (|2^) and ( p6| ) contributes. All critical regions are shown in Fig. |l|, 
which illustrates the analytic structure of both integrands. A given kinematic configuration defines a 
certain compact subset in the depicted C /A, i?/(2A)-plane when the integration variables are varied 
from to 1. In the case of the triangle the curve 

CT„-ti^[B{ti),C{h)] , <ie[0,l] (28) 

represents the integration contour for a given kinematic configuration, i.e. a segment of a parabola. 
In the box case one has 



CBox:{h,t2)^[B{ti,t2),C{ti,t2)] , tje[0,l], (29) 

a family of parabolas. The covering of the corresponding region is multi- valued in general. More 
precisely, B is linear in the integration variables and C quadratic. If thresholds are crossed, C 
typically goes through a minimum and the integration domain is two- valued. If the domains of C^rz 
and Cbox for a given kinematic configuration do not intersect with a boundary of the regions I, II and 
III, the integrand is bounded and no problems arise for a numerical evaluation. For example, in the 
Euclidean region, where all Mandelstam variables are negative, one has always B > Q and C > and 
numerical integration is trivial. An advantageous feature of the sector decomposition (|lj) is that 
it provides bounded integrands in the Euclidean region. For physical processes the domains of Crri 
and Cbox hit the singularities. We display two typical integrands of the triangle graph for a generic 
kinematic configuration in Figs. |^ and |^. In Fig. |^ a situation is shown where only logarithmically 
singular behaviour is present. The imaginary part remains finite in this case. In Fig. ^ one sees 
the square-root singularities when going from region i? < to region I. The real part diverges as 
i? — > from below. The imaginary part diverges as i? ^ from above. This is exactly the behaviour 
expected from eq. ( p7| ) near this boundary. The two logarithmic cusps are located at the boundaries 
of region II. Going from region I to region II leads to a step in the imaginary part. Similar plots are 
obtained by looking at the crossing through regions I and III. In the case of S"^^ one has a family 
of such curves labelled by one of the integration variables. 

A suitable numerical integration routine has to succeed in the critical regions. In principle our 
formulas contain enough information to explicitly detect and classify all possible cases and to make an 
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Figure 1: Analytic regions of the box and triangle integrands. The parabola defines the boundary 
R = 0, the line the boundary A + B + C = 0. Inside regions I, II and III the integrand has an 
imaginary part. The integrable square-root and logarithmic singularities are located at the boundaries 
of these regions as explained in the text. The line segments (a) and (b) stand for the integration 
regions of Figs. 2 and 3, respectively. 
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Figure 2: The integrand of the function S^j.^{6, 4, 1, 1, 1, 1) is plotted for t G [0, 1]. The structure of 
the shown integrand is explained in the text. 



adequate variable transformation at each critical point. Another possibility is to subtract singular 
approximations to the integral, explicitly integrate the singular parts and add the corresponding 
value to the result. In both ways smooth and bounded integral representations can be achieved that 
are free from numerical problems. However, since we have found an iterative numerical method that 
automatically copes with the present integrable singularities, we do not pursue these approaches 
further in this paper. 

We want to close this section by remarking that we implicitly gave a constructive proof of the 
fact that any massive one-loop Feynman diagram has at most integrable singularities of square-root 
or logarithmic type. Further we point out that all the steps of our derivation go through in the 
case of infrared singularities, i.e. in the presence of massless propagators, as the reduction formulas 
are also valid in the massless case and after reduction the IR singularities are isolated in the form 
of triangle functions. The necessary modifications of the given formulas for the remaining infrared 
finite parts are straightforward. 



4 Numerical evaluation 

In the previous sections we showed that all finite scalar A''-point functions can be expressed as 
linear combinations of the "atomic" integrals defined in ( |25| ) and (|26|). In this section we discuss 
two methods that allow to integrate these building blocks numerically without the need for further 
analytic manipulations. 

The characteristics of these elementary integrals were discussed in detail above. In the current 
context we recall first that the integration region is a simple one, namely the unit interval or unit 
square. In this case the optimal approach to achieve high precision rapidly is to use deterministic 
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Figure 3: The integrand of the function S^j,^*{10, 4, 5/2, 1, 1, 1) is plotted for t E [0.05, 0.2] . Integrahle 
square-root and logarithmic singularities are visible as explained in the text. 
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integration rules, particularly for low-dimensional, well-behaved integrals. The case at hand, how- 
ever, is - as also discussed above - complicated by the possible presence of multiple discontinuities 
or integrable singularities. A well suited strategy for these badly-behaved integrands is to divide 
the integration region into subregions, and apply an iterative, adaptive algorithm. This approach 
has been used successfully for irregular high-dimensional integrands in combination with Monte 
Carlo methods p9[|, but has also been explored for lower-dimensional integrands in connection with 
integration rules IpOf . 

Since the integrand of the triangle function (|2^) is only 1-dimensional, we can rely in this case 
on the efficient and highly robust QAGS routine of QUADPACK a widely- used package for 
the numerical computation of definite 1-dimensional integrals. The QAGS algorithm applies an 
integration rule adaptively in subintervals until the error estimate is sufficiently small. The results 
are extrapolated using the epsilon-algorithm, which accelerates the convergence of the integral in 
the presence of discontinuities and integrable singularities. The maximal number of subintervals is 
a fixed input parameter. A maximum of 1000 subintervals proved sufficient to achieve the desired 
relative error of lO"'* in our sample calculations. With this choice, the runtime consumed by triangle 
evaluations turned out to be negligible compared to the one for the box evaluations. 

As succinctly discussed in p2], Section 4.6, evaluating multi-dimensional integrals like function 




(26) poses more difficulties. On the other hand, since the box integral is only 2-dimensional, a 
generalization of the efficient, deterministic method we selected to evaluate the 1-dimensional triangle 
function is suggestive, and several algorithms have indeed been discussed in the literature ]20t . 
Their application, however, would require analytic knowledge of the location of all singularities. 
Hence these algorithms cannot be applied directly to the integral in (|2^)0. We therefore proceed by 
decomposing the potentially singular integrand j{x) into a bounded function hc{x) and a singular 
rest Sc(x) by introducing a cut parameter c > 0: 

/(x) > c 

-c < /(x) < c (30) 
!{x) < -c 

s,(x) fix) - K{x) (31) 

Relying again on integration rules for maximum efficiency, we can then integrate hc{x) with 
DCUHRE 12^, a robust, globally adaptive algorithm applicable to multidimensional, bounded inte- 
grands. To integrate Sc{x) we revert to a well-established, non-deterministic alternative that requires 
no detailed knowledge of the singularity structure of the integrand: Monte Carlo integration ||2^ . 
This advantage, however, is offset by a significantly slower convergence relative to deterministic 
methods. This interplay suggests the existence of an optimal range for the cut parameter c. Be- 
low that range, unnecessarily large, non-singular regions are Monte-Carlo integrated causing slow 
convergence. Above that range, the deterministic routine has to integrate unnecessarily steep peaks 
in the singular regions, necessitating a large number of subdivisions and function evaluations and 
consequently a large workspace due to the global nature of its algorithm. On our systems, the 
workspace limit was about 350MB, allowing for a maximum of 1.5 • 10^ function evaluations for 
DCUHRE. For the calculations shown in the next section we experimented with values for c from 
500 to 50000. The optimal range was 5000 to 10000. The final result is obtained by adding up the 
integrals over hc{x) and Sc{x). As expected, one finds that the obtained results are independent of 
the parameter c. 

We employed an optimised approach to the Monte Carlo integration of Sc{x) that requires a 2- 
dimensional grid covering the integration region. During the integration of bc{x), some evaluations 
of f{x) may return values above or below the cut, i.e. Sc{x) ^ 0. In this case the grid cell that 
contains x is saved. In a second step the integral over Sc(x) is calculated by using all cells with 

^As mentioned above, it is possible to determine analytically the location of the singularities of the box function, 
so in principle, an entirely deterministic integration of the box function is also feasible. 
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Sc{x) 7^ detected in the previous step as seed cells. Sc{x) is integrated in each ceU using crude 
Monte Carlo integration. If a cell result is finite, all neighboring cells are also evaluated. This 
procedure is applied recursively until the region where Sc{x) is non-vanishing is covered. 

Due to its global nature the method described so far requires a potentially large amount of 
memory, and the question arises if a viable "local" alternative with negligible memory requirements 
can be found. To that end, we propose a second, fully recursive approach. Assume the integral Iq 
over a hypercube with volume Vq is to be determined with precision AIq. Starting with volume Vq 
the following procedure is applied recursively: 

1. A value / and error estimate A/ for the integral in the cell of volume V is obtained by applying 

a) an integration rule (ca. 200 integrand evaluations are necessary for a degree 13 integration 
rule ll) 

b) basic Monte Carlo integration with the same number of function evaluations as in a) 

If both results are compatible within errors, the one with the lower error estimate is selected, 
otherwise the result with the larger error is selected. 

2. The tolerable error^in the cell is Almax Aloy^V/Vo. 

3. If A/ < AI„iax no further action is necessary. If AI > Al^ax the cell is divided into n subcells 
of equal volume, and the integrals li in the subcells are determined as in 1. 

1/2 

4. If Aldiv < Al/^/n with Aldiv ■= [X]r=i(^-^0^] ; the procedure is applied recursively to the 
subcells. 

5. If Aldiv > ^I/V^, further subdivision is not advantageous, and / is Monte Carlo sampled 

until AI < Almax- 

This algorithm does not involve a cutoff parameter and was used to double-check the results obtained 
with the first method. For typical box functions with singularities we observe no clear superiority 
of integration rule versus Monte Carlo evaluations, and the added flexibility indeed increases effi- 
ciency. To the best of our knowledge, the combined application of integration rules and Monte Carlo 
sampling has not previously been proposed in the literature. 

Having error estimates for the elementary triangle and box integrals, error estimates for the 
scalar 4-, 5-, and 6-point functions are obtained using standard error propagation. The estimated 
relative error of the results presented in Tables |^ and || is 10"'* or better. 

To conclude, we note that the runtime of our program for the scalar hexagon function depends 
strongly on the complexity imposed by the kinematic configuration and ranges from less than 30 
seconds to many hours. For lower precision the runtime improves considerably. 

5 Results 

To demonstrate the practicality of the approach described above, we calculate the 4-dimensional 
scalar pentagon and hexagon functions for several physical and unphysical kinematic configurations. 
In Table |l| we give numerical results for the pentagon integral. Values are shown for two Euclidean (I 
and II) and two physical kinematic configurations. The latter correspond to scalar integrals occurring 
in the computation of the process 77 — > ttH (pentagon III), and 77 HHH (pentagon IV), see 
Fig. H To work with reahstic mass scales we use Ecms = 800 GeV, mz — 90 GeV, rritop — 175 GeV, 
rriHiggs = 120 GeV. All kinematic invariants are rescaled by {Ecms/'^Y ■ 



^This condition guarantees that the overall error is at most A/q- 
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I 


II 


III 


IV 




-1 


-1 


4 


4 


S23 


-1 


-1 


-1/5 


-7/10 


S34 


-1 


-1 


1/5 


1/10 


S45 


-1 


-1 


3/10 


3/10 


S5I 


-1 


-1 


-1/2 


-1/2 


Si 


-1 


-1 








32 


-1 


-1 








S3 


-1 


-1 


49/256 


9/100 


S4 


-1 


-1 


9/100 


9/100 


S5 


-1 


-5/2 


49/256 


9/100 


mf 


1 


1 


49/256 


49/256 


ml 


1 


1 


49/256 


49/256 


ml 


1 


1 


81/1600 


49/256 


2 


1 


1 


81/1600 


49/256 


2 


1 


1 


49/256 


49/256 


Re 


-0.03542 


-0.03203 


41.33 


3.533 


Im 








-45.96 


-5.956 



Table 1: Four-dimensional scalar pentagon function evaluated for unphysical (1 and II) and physical 
(111 and IV) kinematics. The diagrams defining the kinematics III and IV are depicted in Fig. ^. 
All energies and masses are scaled by Ecmsf^ — 400 GeV. 




(Ill) (IV) 
Figure 4: The diagrams containing the scalar integrals 111 and IV used as benchmarks in Table |l| 
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In Table ^ we display the numerical result for the hexagon integral at the modified|^ symmetric 
point (I) and two physical points (II and III). Point II corresponds to a kinematic situation arising 
in the process 77 — > ttHZ and point III to a kinematic situation in the process 77 HHHH (see 

Fig. D. 
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TT 
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S12 


-1 
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S23 


-I 
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1 /1 n 
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-1/5 


S34 


-i 


1/5 


1 In 

1/5 


S45 


-I 

-i 


/I n 
61 iU 


In 
Z/b 


S56 


-I 

-i 


In 

1/b 


6/ iU 


S61 


-1 


-1/5 


1 / 1 A 

-1/10 


S123 


-1 

-1 


3/11) 


1 / 1 A 

1/10 


S234 


-1 


-1/5 


/ 1 A 

-3/10 


S345 


-0/Z 




U.o»ioyy4d 


Si 


-1 








S2 


1 

- _L 


n 



u 


S3 


-1 


49/256 


9/100 


S4 


-1 


81/1600 


9/100 


S5 


-1 


49/256 


9/100 


S6 


-1 


9/100 


9/100 


9 

ml 


1 


49/256 


49/256 


ml 


1 


49/256 


49/256 


ml 


1 


81/1600 


49/256 


ml 


1 


81/1600 


49/256 


ml 


1 


49/256 


49/256 


ml 


1 


49/256 


49/256 


Re 


0.01353 


-653.8 


-26.93 


Im 





55.56 


48.63 



Table 2: Four-dimensional scalar hexagon function evaluated for unphysical (I) and physical (II and 
III) kinematics. The corresponding diagrams are shown in Fig. |^. All energies and masses are scaled 
by E,^s/2 = 400 GeV. H 





(11) (Ill) 

~ H H~ 
Figure 5: The graphs II and III defining the kinematics for the scalar integrals calculated in Table ||. 



^ The symmetric point for the hexagon function does not obey the nonlinear constraint det G = 0, where G is 
the Gram matrix of the six-point kinematics. The modified symmetric point is designed to satisfy the constraint 
det G = 0, which determines the value for S345. 
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The symmetric and modified symmetric points were used to check our implementation by con- 
firming relations between hexagon and pentagon formulas. To be specific, we considered various 



A^-point functions at the symmetric Euclidean point 


jn=4 
S,sym 




S2 S3 = ~l,ml 


Tn—6 


/r'(s. = 


Si] = -mj = -1) 


rn—4: 
4,sym 




Sy = -m^ = -1) 


Tn=4: 
5,sym 




Sy = = -1) 


jn—A 

b^sym — mod 


■■= ^r'(s5- 


—5/2, else Si — s 


Tn—4 

6,sym-~mod 




= —5/2, else Si — 



-1) = 0.128436 



-1) = -0.0354161 



-0.401140 



-0.0320346 

-1) = 0.013526 



For these special kinematic configurations the reduction simplifies significantly, and one can show 
that the following identities hold: 



A.sym 



jn=& \ 
A,sym } 



Tn=A Tn=i 



(32) 

(33) 
(34) 

Our numerical results - also shown above - fulfill all identities. In addition, we double-checked the 
results using the well-tested program described in Jl7| to calculate multi-loop integrals numerically 
in the Euclidean region. 

The threshold behaviour of the box and hexagon representations is probed in threshold scans, 
which are shown in Figs. ^ and |[ 



Tn=i 

6,sym—7nod 



(Tn=i 
V 5,sym 



^jSym — modJ 



7 




Figure 6: The box and hexagon diagrams which define the kinematics for the scalar integrals used 
in the threshold scans shown in Figures ^ and ^. 

The kinematic configuration for the threshold scan of the box function (Figs. ^ and 0) is 

si2 = {E,„,s/i2mt))\ 323 =- -(mz/(2mt))V2, 
si = S2 = {mz/{2mt)f , S3 = S4 = {nn, / {2mt)f , 
ml = ml = ml = 1/4, ml = {mw / {2mt)Y 

with mf, = 5 GeV and my/ — 80 GeV. 
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The kinematic configuration for the threshold scan of the hexagon function (see Figs. |p and ||) 
is defined by 

Si2 = (£'„„s/(2mf))2, 

S23 = ~S34 = S61 = -Sl23 = -1/10, S45 = -^234 = 2/10, S345 = S345(i?cms), S56 = 3/10, 

Si = S2 = S3 = se 0, S4 = S5 = {mz/{2mt))'^, mf = 1/4. (35) 
Ecms has been varied between 200 and 600 GeV. 




Figure 7: Scan of the 2mt = 350 GeV threshold of the 4-diinensional scalar box function, arising 
from the diagram in Fig. I (a). Details are given in the text. 

We further tested the correctness of the box and pentagon routines by comparing physical results 
with imaginary parts against results obtained using the analytic formulas implemented in ||l^, Q 
and found consistency. 



6 Conclusion 

We have provided explicit representations for the general scalar box, pentagon and hexagon inte- 
grals in terms of n-dimensional triangle and (n-|-2)-dimensional box functions. The advantage of 
such a decomposition is that IR divergences, if present, can be isolated easily, such that an efficient 
book-keeping of IR poles is achieved. The remaining integrals can be evaluated by setting e = 0. In 
our approach, the finite triangles and 6-dimensional box functions are represented as one- respec- 
tively two-dimensional parameter integrals in a form which is convenient for numerical integration. 
Although in principle all these integrals can be evaluated analytically by applying scalar reduction 
formulas, it is clear that the expression for the general hexagon function would contain a huge num- 
ber of dilogarithms which would give rise to nontrivial cancellations. Representations with a smaller 
number of "atoms" seem to be preferable from this point of view. 

This motivated us in the present work to follow a numerical approach at an earlier stage of 
the calculation, where the "atoms" are the n-dimensional triangle and the (n-|-2)-dimensional box 
function instead of logarithms and dilogarithms. Focusing on the massive case we studied in detail 
the singularity structure of our one (two)-dimensional integral representations for the triangle (box) 
functions. Real and imaginary parts are explicitly separated in this approach which avoids to deal 
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80 



50 




200 300 400 500 600 200 300 400 500 



Ecms [GeV] E,^, [GeV] 

Figure 8: Scan of the 2mt = 350 GeV threshold of the 4-diniensional scalar hexagon function, see 
Fig. ^D. The kinematics is defined in the text. 

with infinitesimal quantities. We have shown for physical kinematics that the integral representations 
for triangle and 6-dimensional box have an analogous singularity structure which can be treated 
transparently for both in the same way. 

In principle our formulation contains enough information to produce bounded and smooth in- 
tegrands by adequate variable transformations or subtractions at the critical points. However, it 
seemed an interesting question to us whether such somewhat cumbersome manipulations can be 
avoided by directly evaluating the integral representations. To our best knowledge no general nu- 
merical algorithm has been provided in the literature so far for the relatively complicated singularity 
structure present in the 2-parameter integral representation of the 6-dimensional box. By combining 
deterministic integration methods, adequate for the smooth part of the integrand, and Monte Carlo 
techniques for the singular regions in an iterative way, we succeeded in numerically evaluating the 
integrals directly. Two different numerical integration methods have been presented. The correct- 
ness and efficiency of our integration strategies has been demonstrated by comparing our results to 
known results if available. For the hexagon integral, which could only be tested indirectly, we have 
performed several consistency checks. We have shown that a stable result can be obtained with our 
method when internal thresholds are probed. 

A natural question to ask next is if a fully numerical approach to the evaluation of multi-particle 
production at one-loop is feasible. Analytic calculations are generally hampered by the enormous 
complexity generated when reducing integrals with nontrivial numerators to scalar integrals. It is 
conceivable that the critical reduction steps could be avoided to a large extent if all or some groups 
of Feynman diagrams are treated numerically, such that one is dispensed from doing the full tensor 
reduction. We note in this respect that the numerical difficulties stem entirely from the denominators 
of the integrals, whereas tensor integrals only introduce additional parameters in the numerators, 
which pose no problem. 

Following our strategy one is always able to separate real and imaginary parts, after having done 
the trivial integrations analytically. The remaining integrals will always have at most singularities 
of square root and/or logarithmic type, which can be integrated directly with our methods. Thus 
our findings can be viewed as a step towards a complete numerical approach to calculate multi-scale 
processes at one loop. 
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